function idat = fsb_interpolate_filter_scans(idat,sandbox)

% FSB : interpolate out-of-trial periods in concatenated data and filter data
%
% EXAMPLE:
% idat = fsb_interpolate_filter_scans(idat,sandbox)
%
% INPUT:
% idat: 4-D image data
% sandbox: sandbox experiment struct
%
% OUTPUT:
% idat: interpolated 4-D image data
%
% CALLED BY:
% FSB.m
%
% NOTES:
% It is recommended to use this function instead of only the interpolation
% function
% Function accepts concatenated data and outputs concatenated data, too
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $Revision 1.0
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Check for correct dataset size
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

d = size(idat);
if size(d)<4;
    disp('No 4D dataset, aborting');
    return
end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Run function with sandbox input
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if min(sandbox.intrial(:,1))<1
    disp('please concatenate trials first!')
    return
end

if size(sandbox.intrial,2)<4
    sandbox.intrial(:,4)=1;
end

for scannum = 1:max(sandbox.intrial(:,4));

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Get the respective segment of image data and inflate it
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    scanindex = find(sandbox.intrial(:,4) ==scannum);
    idat_seg = single(idat(:,:,:,scanindex));
    idat_dummy = zeros(d(1), d(2), d(3), max(sandbox.intrial(scanindex,7)));
    subscanindex = sandbox.intrial(scanindex,7);
    idat_dummy(:,:,:,subscanindex) = idat_seg;
    
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Get the respective segment of experiment data and inflate it
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    intrial_dummy = zeros(max(sandbox.intrial(scanindex,7)),7);
    intrial_seg = sandbox.intrial(scanindex,:);
    intrial_dummy(subscanindex,:) = intrial_seg;
    intrial_dummy(:,4)=scannum;   
    intrial_dummy(:,5)= scannum;
    intrial_dummy(:,6) = scannum;
    intrial_dummy(intrial_dummy(:,2)==0,2) = 1;

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Get minimum and maximum of scans
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    trials = unique(intrial_dummy(:,3));
    trials(trials ==0) = [];
    trialnum = length(trials);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Get overall trial indices
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    trialind = find(intrial_dummy(:,2)>0 & intrial_dummy(:,1)>0); % This is the source of the problem, presumably
    trialind2 = intrial_dummy(trialind,7);
    dummy2 = zeros(max(intrial_dummy(:,7)),1);
    dummy2(trialind2) = 1;
    dummytest2 = bwlabeln(dummy2);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Calculate average length of trials
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    for volnum = 1:trialnum;
        trialind3 = find(dummytest2==volnum);
        testi(volnum) = length(trialind3);
    end

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Find indices for last volumes of each period
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    
    endindex = [];
    startindex = [];
    
    trialmax = max(intrial_dummy(:,3));
    trialminindex = find(intrial_dummy(:,3)>0);
    trialmin = min(intrial_dummy(trialminindex,3)>0);

    for trial_num = trialmin:trialmax;
        temptest = find(intrial_dummy(:,3)==trial_num);
        maxvol = max(intrial_dummy(temptest,7));
        endindex = [endindex maxvol-1 maxvol];
        minvol = min(intrial_dummy(temptest,7));
        startindex = [startindex minvol minvol+1];
    end

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Create indices of trial beginning and end
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    trialindex1 = find(intrial_dummy(:,2)<3 & intrial_dummy(:,1)~=0);
    % Startindex
    trialindex1 = intrial_dummy((trialindex1),7);
    try
        if trialindex1(:)~=startindex(:)
            disp('Mismatch');
        end
    catch
        disp('')
    end

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Create dummy arrays to place data into and label them
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    trialindexdummy1 = zeros(max(intrial_dummy(:,7)),1);
    trialindexdummy2 = trialindexdummy1;
    trialindexdummy1(startindex) = 1;
    trialindexdummy2(endindex) = 1;

    test1 = bwlabeln(trialindexdummy1);
    test2 = bwlabeln(trialindexdummy2);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Go through trials and figure out what needs to be interpolated
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    g = waitbar(0);

    volind1 = {};

    for volnum = 1:trialnum;

        % Graphics output
        detstring = ['Calculation artefact interpolation for scan' num2str(scannum) ', trial ' num2str(volnum) ' of ' num2str(trialnum)];
        disp(detstring)
        waitbar(volnum/trialnum,g,detstring);
        
        % Actual calculations
        volind1{volnum} = find(test1==volnum);
        startvol = idat_dummy(:,:,:,volind1{volnum});
        startvolt{volnum} = mean(startvol,4);
        volind2 = find(test2==volnum);
        
        % Control condition
        try
        endvol = idat_dummy(:,:,:,volind2);
        catch
            disp(volind2);
        end
        endvolt{volnum} = mean(endvol,4);

    end
    close(g);

    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    % Interpolate Artefact volumes
    %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

    g = waitbar(0);
    if trialnum>1
        for volnum = 2:trialnum;

            % Graphics and feedback output
            detstring = ['Interpolating artefacts for scan' num2str(scannum) ', trial ' num2str(volnum) ' of ' num2str(trialnum)];
            disp(detstring)
            waitbar(volnum/trialnum,g,detstring);
            
            %  determine volume difference
            voldiff = startvolt{volnum} - endvolt{volnum-1};
            volind2 = find(test2==volnum-1);

            % Control condition
            try
                artlength = volind1{volnum}(1)-max(volind2);
            catch
                disp('ATTENTION: Error while determining length of artefact period')
                disp(volind2);
            end

            artsec = repmat(endvolt{volnum-1},[1 1 1 artlength-1]);
            try
            voldiff2 = voldiff/artlength;
            catch
                disp('')
            end
            
            voldiff2 = repmat(voldiff2,[1 1 1 artlength-1]);
            artgrad2 = 1:artlength-1;
            artgrad2 = repmat(artgrad2,[size(voldiff2,1) 1 size(voldiff2,2) size(voldiff2,3)]);
            artgrad2 = permute(artgrad2,[1 3 4 2]);
            artgrad3 = artgrad2.*voldiff2;
            artgrad3 = artgrad3+artsec;

            try
                idat_dummy(:,:,:,max(volind2)+1:volind1{volnum}(1)-1) = artgrad3;
            catch
                disp('ATTENTION: Error while interpolating volumes')
            end

        end

        try
            endind = find(test2==volnum);
            endseg = max(endind)+1:size(idat_dummy,4);
            idat_dummy(:,:,:,endseg) = repmat(endvolt{volnum},[1 1 1 length(endseg)]);
            startseg = 1:min(volind1{1})-1;
            idat_dummy(:,:,:,startseg) = repmat(startvolt{1},[1 1 1 length(startseg)]);
        catch
            disp('Issues with start or end segment')
            endind = [];
        end

        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        % Filter data with selected filter
        %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

        %sandbox.filttype = 'Stop';
        idat_seg = fsb_filter_temporal_single(idat_dummy,sandbox);        
        idat(:,:,:,scanindex) = idat_seg(:,:,:,trialind2);

    else
        volnum =1;
        idat(:,:,:,scanindex) = idat_dummy(:,:,:,trialind2);
        disp('Not filtering data, only one trial in scan');
    end

    disp(['...done with scan ' num2str(scannum)]);
    close(g);
end

end
